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^ . ABSTRACT 

The population of 7-ray emitting millisecond pulsars (MSPs) is studied by 

2 ' using Monte-Carlo techniques. We simulate the Galactic distributions of the 

+3 ; 

^ ' MSPs, and apply the outer gap model for the 7-ray emission from each sim- 

! 

ulated MSP. We take into account the dependence of the observed 7-ray flux 
on the viewing angle and inclination angle, which is the angle between the 
' rotation axis and the magnetic axis, respectively. Using the sensitivity of the 

' six-month long observation of the Fermi telescope and radio sensitivities of 

(N ■ 

I are detected within our simulation. The statistical properties of the simulated 



existing pulsar surveys, 9-13 radio-selected and 22-35 7-ray-selected pulsars 



population are consistent with the Fermi observations. Scaling the observed 

sensitivity cx a/T, where T is the length of observation time, the present model 
• 1—1 

I predicts that over the 5-year mission Fermi would detect 15-22 radio-selected 

H ■ 

. 7-ray MSPs, and 95-152 7-ray-selected MSPs. Our simulation also predicts 

that about 100 (or 200-300) 7-ray MSPs with a flux larger F ^ 10"" erg/cm^s 
(or 5 X 10"^^ erg/cm^s) irradiate the Earth. With the present sensitivities of 
the radio surveys, most of them are categorized as 7-ray-selectcd pulsars, in- 
dicating that most of the 7-ray MSPs have been missed by the present Fermi 
observations. We argue that the Galactic Fermi unidentified sources located 
at high latitudes should be dominated by MSPs, whereas the sources in the 
galactic plane are dominated by radio-quiet canonical pulsars. We want to 
emphasize that the predicted number of radio-loud 7-ray MSPs depends on 
the sensitivities of radio surveys and that is can be increased, for example, 
from 15-22 to 26-37 if the radio sensitivity is improved by a factor of 2. 
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1 INTRODUCTION 

The millisecond pulsars (MSPs), which have a rotation period P ~ 0.002 — 0.1 s and a 
stellar magnetic field ~ 10^^^ G, are classified as a different generation from the canonical 
pulsars, which have P ~ 0.03 — 10 s and Bg ~ 10^^^^'^ G. The Fermi LAT has detected 7-ray 
emissions from about 60 pulsars in just two year observation, including 9 millisecond pulsars 
(Abdo et al. 2010a,b,c, 2009a,b; Saz Parkinson et al. 2010). Furthermore, the detection of 
radio MSPs from about 20 unidentified Fermi point sources (e.g. Ray 2010; Caraveo 2010; 
Ranson et al. 2011) has been reported. The MSPs are now recognized as one of the major 
populations of the Galactic 7-ray sources. 

The particle acceleration and high-energy 7-ray radiation processes in the pulsar mag- 
netosphere have been studied with the polar cap model (Ruderman & Sutherland 1975; 
Daugherty & Harding 1982, 1996), the slot gap model (Arons 1983; Muslimov & Harding 
2003; Harding et al. 2008) and the outer gap model (Cheng, Ho & Ruderman 1986a,b; 
Hirotani 2008; Takata, Wang & Cheng 2010a) respectively. All models have assumed that 
the charged particles are accelerated by the electric field along the magnetic field lines, and 
that the accelerating electric field arises in the charge deficit region, where the local charge 
density deviates from the Goldreich- Julian charge density (Goldreich & Julian 1969). The 
polar cap model assumes that the acceleration region is near the stellar surface. On the 
other hand, the outer gap model and the slot gap model assume a strong acceleration region 
extending to the outer magnetosphere. For the canonical 7-ray pulsars, several access of 
observational evidence have been proposed for the outer magnetospheric origin of the 7-ray 
emissions (Aliu et al. 2008; Abdo et al., 2009c). Because the strength of the magnetic field at 
the outer magnetosphere of the MSPs is similar to that of canonical pulsars, B ~ 10^"^ G, it 
is expected that if the 7-ray emissions of the MSPs originates from the outer magnetosphere, 
the emission characteristics (e.g. pulse profile and spectrum) will be similar to those of the 
canonical pulsars. 

Romani & Watters (2010) and Watters & Romani (2011) studied the pulse profiles of 
the canonical pulsars observed by Fermi. They computed the pulse profiles predicted by the 
outer gap model and slot gap model, and argued statistically that the outer gap geometry 
is more consistent with the observations than the slot gap geometry. Venter, Harding & 
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Guillemot (2009) fit the pulse profiles of the Fermi detected MSPs with the geometries 
predicted by the different high-energy emission models. They found that the pulse profiles 
of two out of eight millisecond pulsars cannot be fitted by the geometries of either the outer 
gap or the slot gap models. They proposed a pair-starved polar cap model, in which the 
multiplicity of the pairs is not high enough to completely screen the electric field above the 
polar cap, and the particles are continuously accelerated up to high altitude over the entire 
open field line region. Thus we see that the 7-ray emission mechanism and the emission site 
in the MSP's magnetosphere have not been satisfactorily explained up to now. 

The increase of population of detected 7-ray emitting pulsars allows us to perform a 
more detailed population study of the high-energy emissions from the pulsars. A comparison 
between the simulated and observed distributions of the various pulsar characteristics (e.g. 
rotation period and 7-ray flux) will be useful to test the high-energy emission model. With 
the canonical 7-ray pulsars, Takata, Wang and Cheng (2011) have studied the population 
predicted by the outer gap model. They predicted that with the sensitivity of the six- 
month long Fermi observation, about 100 7-ray emitting canonical pulsars can be detected, 
suggesting that the present observations have missed many 7-ray emitting pulsars. For the 
MSPs, Story, Gonthier and Harding (2007) have studied the population of 7-ray MSPs 
with the slot gap accelerator models and predicted the Fermi observations. They predicted 
Fermi will detect 12 radio-loud and 33-40 radio-quiet 7-ray MSPs. 

The population study of the MSPs is also important for an understanding the unidentified 
7-ray sources detected by Fermi. In the Fermi first source catalog, there are several hundred 
unidentified steady point sources (Abdo et al. 2010b). Takata et al. (2011) have argued that 
the Fermi unidentified sources located in the high Galactic latitudes cannot be explained 
by the galactic distribution of the canonical pulsars. The millisecond pulsars are possible 
candidates for the high Galactic Fermi unidentified sources. 

In this paper, we develop a Monte-Carlo study of the population of the 7-ray emitting 
MSPs predicted by the outer gap model. In particular, we only study the population of MSPs 
in the Galactic field (not in the globular clusters). Following the previous studies (e.g. Story 
et al. 2007; Takata et al. 2011), we will perform a Monte-Carlo simulation for the Galactic 
population of the MSPs (section [21). In section [231 we will discuss our 7-ray emission model, 
including dependence of the 7-ray flux on the viewing geometry. In section [3l we present the 
results of our simulation. In particular, we will compare the simulated population with the 
Fermi six-month long observation (section 13.21) . In section [331 'we discuss the possibility of 
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MSPs as the origin of the Galactic Fermi unidentified sources. We will discuss the results 
of our simulation in section HI 



2 MONTE-CARLO SIMULATION 
2.1 Galactic population 

In this study, we assume that all MSPs are born through the so called recycled process in 

the low-mass binary systems, with a birth rate of 10"^ ~ 10~^ yr"'^ (Lorimer et al. 2005; 

Lorimer 2008). The birth location is determined by the spatial distributions proposed by 

Paczynski (1990), 

aRe-^/^-pi? 
Pb{R) = , 

Pz{z) = ^e-l^l/^-^ (1) 

^cxp 

where R is the axial distance from the axis through the Galactic centre perpendicular to 
the Galactic disk and z is the distance from the Galactic disk. In addition, -Rexp = 4.5 kpc, 
an = [1 - e--f^™''-/-f^'=-f (1 + Rmax/ Rexp)]~'^ with R^ax = 20 kpc and Zexp = 200 pc (Story et 
al. 2007). 

To obtain the current position of each simulated MSP, we solve the equations of motion 
from its birth to the current time. The equations of motion are given by 
dR^ vl d^tnt ^2) 



df^ R OR ' 
dz^ _ d^tot 
dt^ ~ dz ' 

and 

Rvtj, = constant. (4) 

Here is the azimuthal component of the velocity, ^tot = ^sph + ^dis + is the total 
gravitational potential, where ^sph, ^dis and $/j are the spheroidal, the disk and the halo 
components of the Galactic gravitational potential, and are given by 

MR, ^) = - I (5) 

where i = sph and dis. For the spheroidal component, asph = 0, bgph = 0.277 kpc and 
Mgph = 1-12 X 10^'^Mq. For the disk component, adis = 3.7 kpc, bais = 0.20 kpc, and 
Mdis = 8.07 X 10^^ Mq, while for the halo component 
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CM 



(6) 



where = 6.0 kpc Mc = 5.0 x 10 Mq, respectively (c.f. Burton & Gordon 1978; Binney & 
Tremaine 1987; Paczynski 1990). The Lagrangian in units of energy per unit mass is given 
by 

L = ^^^-t„(fl..). (7) 

where v is the velocity of the pulsar. 

For the distribution of the initial velocities of the MSPs, we assume a Maxwellian distri- 
bution with a characteristic width of = 70 km/s (c.f. Hobbs et al. 2005), namely, 



Mv) = Ji,-,e-^ (8) 

For the azimuthal component, the circular velocity due to the Galactic gravitational potential 
field at the birth position of the MSPs is taken into account, and it is calculated from 



R 



OR OR OR 



(9) 



2.2 Pulsar characteristics 



It is widely accepted that the MSPs are so called recycled pulsars, which were spun up by 
accretion of the matter from the low mass companion star. The rotation period of the newly 
born MSP is related to the history of the accretion onto the neutron star. In particular, 
its rotation period in the accretion stage may be related to the equilibrium spin period. 
Peg ~ l.TB^^"^ Mi^^"^ R^Q^"^ Mi^^"^ , where is the neutron star magnetic field in units of 
10*^ Gauss, Mi5 is the accretion rate in units of 10^^ g/s, R^ is the neutron star radius in 
units of 10® cm and Mi. 4 is the neutron star mass in units of 1.4 solar mass (Frank, King & 
Raine, 2002). However, the description of the transition from an accretion powered to the 
rotation powered phase is not well understood due to the complexities in the description 
of the interaction between the magnetosphere of a neutron star and its accretion disk. 
Furthermore, the transition, for which a rapidly decreasing accretion rate is required (e.g. 
Jeffrey 1986), may be facilitated by different processes between the long and the short orbital 
period systems. For the long orbital period systems, a red giant companion can detach from 
its Roche lobe as its envelope is exhausted. On the other hand, for the short orbital period 
systems, a possible mechanism leading to the sudden decreas in the accretion flow is the 
operation so called "propeller effect" (Campana et al. 1998), or dissolution of the disk by 
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7-ray irradiation in the quiescent stage (Takata, Cheng & Taam, 2010b). The initial period 
of the MSPs will depend on the characteristic of the orbital motion and/or the termination 
mechanism of the accretion flow. 

In the present study, we assume that all MSPs are not directly produced by supernova 
explosions, but that they are born by the recycled process in a low-mass binary system. On 
this assumption, the true age of the MSPs, defined by the time since the neutron star was 
born in the supernova explosion, is different from the spin down age (i.e. P/2P), and from 
the age defined by activating the rotation powered activity. Because the time scale of the 
recycled process, including the decay of the neutron star's magnetic field and the accretion 
process, is not understood well, it is very difficult to calculate the present period and the 
magnetic field strength of MSPs from their initial distributions and the true age of MSPs. 

Allowing for uncertainties, we use the observed distribution of the radio MSPs to assign 
"current" pulsar characteristics (e.g. rotation period, magnetic field) for each simulated MSP, 
instead of modeling the initial distribution. In the usual Monte-Carlo studies for the Galactic 
population of pulsars (e.g. Story et al. 2007; Takata et al. 2011), the initial period is assigned 
for each simulated pulsar and then the current position and rotation period are obtained. In 
the present study, on the other hand, we (1) generate the simulated MSPs with a constant 
birth rate over 10 Gyr, (2) obtain the current position as described in section 12.11 and (3) 
assign the pulsar parameters following the observed distributions. Specifically, we assign 
the period time derivative (P) and the stellar magnetic field (Pg) for each simulated MSP 
following the observed P — Bg distribution, where we use the intrinsic value after removing 
the Shklovskii effect (see below). In this paper, we denote as P^ the strength of the magnetic 
field at the magnetic equator, namely, Pg = 3.2 x 10^^ vPP. The current rotation period 
and the spin down age of simulated MSPs are calculated from 

P_3 = 0.97p2pr2'o ms (10) 
and 

r = 1.5 X 109p_3P2o' yr (H) 
respectively (Lyne & Graham-Smith, 2006). Here P_3 and P--20 are the rotation period in 
units of 0.001 s and the period time derivative in units of 10~^°, respectively. The distribu- 
tions, which are used in the present simulation, are represented for the various characteristics 
of the simulated MSPs in Figure [1] by the dashed lines. Figure [1] also presents the observed 
distributions (shaded-histograms, Manchester et al. 2005). 
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Shklovskii (1970) argued that the Doppler shift resuhing from the transverse motion of 
a pulsar makes a positive contribution to the pulsar's period time derivative as (Manchester 
1999), 

p. - ^, (12) 

where fi is the proper motion of the MSP. For the MSPs whose period time derivative is 
extremely small, the Shklovskii effect may significantly increase the period time derivative. 
For canonical pulsars, on the other hand, the period time derivative is large, e.g. P ~ 
10"^"^ — 10^^^, which is much greater than that caused by the Shklovskii effect. 

We note that the present procedure can be applied if the Galactic distribution of the 
MSPs does not depend on the age of the rotation powered MSPs (and spin down age); that 
is, the Galactic MSP's populations of 1 Gyr and of 10 Gyr MSPs, for example, are described 
by the same distribution. With the present assumption that all MSPs are born through 
the recycled process, which has a characteristic time scale completely independent of the 
spin down age, the Galactic distribution is independent of the distribution of the spin down 
age. With a typical velocity of the observed MSPs , f ~ 70 km/s, it is expected that the 
typical displacement of MSPs with an age, r ^ 100 Myr, becomes larger than the size of 
the Galaxy. However, with the slow velocity, the MSPs remain bound to the Galaxy and 
hence their Galactic distribution does not depend on the age of the rotation powered MSPs. 
On the other hand, for canonical pulsars (in particular for the 7-ray pulsars), the Galactic 
distribution depends on the true age, which is almost equal to the spin-down age. This allow 
us to calculate the present period distribution from the initial distribution. 

2.3 Radio emissions 

Using the empirical relations among the radio luminosity, rotation period, and period time 
derivative, the distribution of the radio luminosity at 400 MHz is expressed by (Narayan & 
Ostriker 1990) 

PU,, = 0.5X'e\ (13) 

where A = 3.6[logio(L4oo/ < L400 >) + 1.8] with L400 = r]10^-^'^P^^^/P^, and L400 is the 
luminosity in units of mJy kpc^. Here 77 is a scaling factor to adjust the observed distribution, 
and ?7 = 1 for the canonical pulsars. In the present simulation for the MSPs, we find that 
rj ~ 0.05 can explain the distribution of the observed radio luminosity of the MSPs. The 
radio flux on Earth is given by 6*400 = L^Qo/d"^, where d is the distance to the MSP. We scale 
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the simulated 400 MHz luminosity to the observational frequency using a typical photon 
index ~1.8 (Kramer et al. 1998). 

We also take into account the beaming effect of the radio emissions. The half-angle, 
which is measured from the magnetic axis, of the radio emission cone of the MSPs does not 
depend on the frequency, and is approximately described by (Kramer & Xilouris, 2000), 



This emission can be detected by observers with a viewing angle between max(0°,a — u) 
and min(a + w, 90°), where a is the inclination angle between the rotation axis and the 
magnetic axis. 

We use the ten radio surveys (Molongo 2, Green Band 2 and 3, Arecibo 2 and 3, Parks 
1, 2 and MB, Jordell Bank 2 and Swinburne IL), whose system characteristics are listed 
in Table 1 of Takata et al. (2011) and the references therein. To calculate the dispersion 
measure, we use the Galactic distribution of electrons obtained by Cordes & Lazio (2002). 

2.4 7-ray emission model 

2.4-1 Observed 'y-ray flux 

In this paper we apply the outer gap accelerator model (Cheng, Ho & Ruderman 1986a,b; 
Zhang & Cheng 2003; Takata et al. 2010a) for the millisecond pulsars. In the outer gap, the 
electrons and/or positrons are accelerated up to a Lorentz factor of F ^ 10^ by the electric 
field along the magnetic field line. These accelerated particles can emit 7-ray photons of 
several GeV through the curvature radiation process. Assuming the force balance between 
the electric force and the radiation drag force, the Lorentz factor is given by 



where Rc is the curvature radius of the magnetic field lines and E\\ is the accelerating electric 
field. The accelerating electric field in the outer gap is approximately given by (Cheng et al. 
1986a,b; Cheng, Ruderman & Zhang 2000) 



where Ric = cP/2ii is the light cylinder radius and /(r) is the gap thickness divided by 
the light cylinder radius {Ric = cP/2'k) in the poloidal plane. In addition, z is the height 
measured from the last-open field line in the poloidal plane, and z = Zhis the upper boundary 
of the gap. We note B{r)f{rY ~constant along the field line for the dipole field geometry. 



(14) 




(15) 




(16) 
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The spectrum of the curvature radiation emitted by an individual particle may be written 

as 

^c(i^7,0 = ^i^(^), (17) 

where x = E^/Ec with E^ = ShcT^/AirRc and 

POO 

E{X) = X / Kr,/3{t)dt, 

where K^j^ is the modified Bessel function of order 5/3. If the 7-ray beam points toward an 
observer, the observer will measure a phase-averaged spectrum of (Hirotani 2008), 

dE.^T^J Nin)Rc{n)Pc{E„ r^)dA, (18) 

where N = B/Pce is the particle number density, represents the radius from which the 
emission can be measured by the observer, dA is the cross section of the gap perpendicular 
to the magnetic field lines. We assume that the outer gap extends in the azimuthal direction 
with 6(f) ^ 71 radian. The integrated energy fiux between 100 MeV and 300GeV can be 
calculated from 

r300GeV 

-^7,>100 MeV = / -j-=;-dEy. (19) 

JlOOMeV dE^ 

One can show that the 7-ray fiux described by equation (fT9|) approximately satisfies 
-^7,>ioo MeV oc /{RicYLgd, where Lgd = 4(27r)^i?^i?^/6c^P^ is the spin down power. If the 
typical energy of a curvature photon E^ satisfies 100 MeV ^ E^. ^ 100 GeV, the power 
radiated by a single charged particle is / P^dE^ ~ 2e^cr^/3i?^. Furthermore, if we estimate 
the total number of charged particles in the gap in terms of values near the light cylinder as 
^aap ~ ^^rRic^MRic), we obtain 



S,>100 MeV ~ -^^^^aap ^ « ^2p4g3 °^ ^ ^^-^/^ ' ^^^^ 



where we have used T ~ (3i^^2^E||/2e)l/^ E^ ~ B{Ric)P{Ric) and 5A ~ f5<pRl. The 7-ray 
fiux is proportional to the cube of the fractional gap thickness and the spin down power. 

We also assume as a zeroth order approximation that the gap current is of the order of the 
Goldreich- Julian value over the full width. The detailed calculation of the outer gap model 
(e.g. Takata & Chang 2007) gives the current distributions in the direction perpendicular 
to the magnetic field lines. On the other hand, by fitting the observed gamma-Taj spectra 
of all mature pulsars detected by Fermz, Wang, Takata & Cheng (2010) find that the total 
current in the gap is of the order of the Goldreich- Julian value. We expect that as long as the 



10 J. Takata, Y. Wang and K.S. Cheng 

total power of the 7-ray radiation from the outer gap accelerator is concerned, the uniform 
current distribution with the Goldreich- Julian value would not be a bad approximation. 



2.4-2 Thickness of the outer gap 

Let's discuss now the thickness of the outer gap accelerator. Zhang & Cheng (2003) estimated 
the gap thickness for the millisecond pulsars when the photon-photon pair-creation process 
controls the gap activities. They have argued that the gap thickness is controlled by the 
photon-photon pair-creation process between the 7-rays emitted in the outer gap and the 
X-rays from the stellar surface, where the multiple magnetic field dominates the global dipole 
field (Ruderman 1991, Arons 1993). They estimated the typical gap thickness divided by 
the light cylinder radius at the light cylinder as 

fzc = = 7.0 X 10-^ PT'B^'^'Srr. (21) 

J^lc 

Here 6rr, is the distance (in units of 10^ cm) from the stellar surface to the position where the 
local magnetic field is equal to the dipole field. In the following we assume that Sr^ ~ 1 — 10. 

Takata et al. (2010a) proposed the outer gap model, in which the gap thickness is de- 
termined by the magnetic pair-creation process near the stellar surface. They argued that 
the returning particles, which were accelerated in the gap, will emit photons with an energy 
meC^/af ~ lOOMeV by curvature radiation near the stellar surface, where mgC^ and a/ are 
the electron rest mass energy and the fine structure constant, respectively. The lOOMeV 
photons form pairs by the magnetic pair creation process. The photon multiplicity is easily 
over 10^ for each incoming particle. For a simple dipole field structure, all pairs should move 
inward and cannot affect the outer gap accelerator. However Takata et al. (2010a) have ar- 
gued that if the magnetic field lines near the surface, instead of being nearly perpendicular 
to the surface, are bending side-ways due to the strong local field, the pairs created on these 
local magnetic field lines can have an angle bigger than 90°, which results in an outgoing 
flow of pairs. With this model, the estimated fractional gap thickness is given by 

fm = ~ = 0.8A1P/1 s)V^ (22) 

Kic Kp 

where Rp ~ Rs{Rs/ RicY^'^ is the polar cap radius, K ~ B^i2Sj is the parameter character- 
izing the local magnetic field properties, 5m,i2 and S7 are the local magnetic field in units of 
lO^^G and the local curvature radius in units of lO^cm, respectively. By fitting the radiation 
properties of the 7-ray pulsars observed by Fermi, they find that K 2 for the canonical 
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pulsars, while ~ 15 for the millisecond pulsars. When the fractional gap thickness fm 
is smaller (or larger) than /^c, the magnetic pair-creation (or photon-photon pair-creation) 
process controls the gap thickness. 

Finally we discuss the maximum fractional thickness, fmaxi for the active outer gap ac- 
celerator. Zhang & Chang (1997) have argued that a pulsar with the fractional gap thickness 
larger than unity, fzc > 1, is not active, because the pairs are not created in the gap by 
the photon-photon pair-creation process. However, the outer gap accelerator can exist only 
between the last-open field lines and the critical field lines, which are defined by those mag- 
netic field lines that have the null point of the Goldreich- Julian charge density at the light 
cylinder. Therefore, we may define the maximum gap thickness as 

fcrit = 2) ' (^^) 

where 9p and 9c are the polar angles of the last-closed field line and of the critical field line 
on the stellar surface, respectively. For the pure dipole field, we obtain 



= a -|- sin ^ 



and 

6c = a + sin~^ 



sm{9ic - a) — ^sin6',, 



ic 



sin(6'„ — a) [ sin 9n 



R. 



Ic 



(24) 



(25) 



respectively, where 

1 /-3 - \/9 + 8tan2a^ 



4 tan a 



9ic = tan 
and 

_i / 3 tan a + \/8 + 9 tan^ a \ 

9n = tan . 

V 2 J 

In this paper, we examine two extreme cases corresponding to /max=l and fmax = fcrit 
respectively. The results for these two cases may give a range of uncertainty of the present 
theoretical predictions. 

2.4-3 Dependence on the inclination and viewing angles 

The observed characteristics of the 7-ray emissions depend on the viewing angle and the 
inclination angle between the rotation axis and the magnetic axis. With the geometry of 
the outer gap model, stronger emissions are emitted toward a viewing angle ^ ~ 90°, and 
the observed flux tends to decrease as the viewing angle closes to the rotation axis, where 
^ = 0° or ^ = 180°. To take into account this effect, we calculate the 7-ray flux, described 
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in section 12.4. ![ as a function of the viewing angle and of the inchnation angle. In the 
calculation, we apply the vacuum dipole field geometry and assumed that the outer gap 
extends from the null charge point of the Goldreich- Julian charge density to r = l.'bRic or 
p = Ric, where p is the axial distance from the rotation axis. We also assume for simplicity 
no azimuthal dependence of the gap structure. 

Figure [2] shows the dependence of the 7-ray flux as a function of the viewing angle (C,) and 
of the inclination angle [a). The vertical axis represents the fractional 7-ray flux, which is 
defined as the 7-ray flux divided by f^Lsd/d"^- The different lines correspond to the different 
inclination angles. The results are obtained for a rotation period P = 3 ms and the stellar 
magnetic field Bg = 3 x 10^ Gauss. We can see in Figure [2] that the calculated flux tends to 
decrease as the viewing angle and the inclination angles decreases from 90°. 

We can see that the general trend of the relation between the fractional 7-ray flux and 
the viewing angle and the inclination angle, seen in Figure 121 is maintained for different sets 
of the rotation periods and magnetic fields. In this study, for simplicity, we use the relations 
in Figure [2] for all sets of the rotation period and of the magnetic field. We assume that the 
inclination angle and the viewing angle are randomly distributed. 

However, we do not take into account the dependence of the spin down power on the 
inclination angle. The pulsar spin-down can result from both the current braking torque and 
the magnetic dipole radiation. The former can spin-down the pulsar even if the inclination 
angle is zero. According to the analysis of the force-free magnetosphere done by Spitkovsky 
(2006), the spin down power depends on the inclination angles as Lg^ oc 1 -|- sin^a; in 
other words, Lg^ changes only by a factor of 2 with the inclination angle. With this small 
variation, we may expect that the general properties of the simulated distributions for the 
MSP's characteristics (such as the period and the magnetic field) are not affected much by 
dependence on the inclination angle. 

3 RESULTS AND DISCUSSIONS 

In the Monte-Carlo simulations, we generated 3 x 10^ MSPs with a constant birth rate 
over 10 Gyr. About 2100 radio MSPs were detected by the simulated radio surveys. Scaling 
the simulated population to the number of radio MSPs in the ATNF catalog, ~ 80, we 
obtain ~ 10~^ per year as the predicted averaged birth rate, which is consistent, within the 
uncertainties, with the other population studies (e.g. Ferrario & Wickramasinghe 2007). 
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3.1 Galactic distribution of MSPs 

In this section, we discuss the consistency of the observed and simulated distributions of 
the radio MSPs. First, Figures [3] (a) and (b) compare the observed (shaded histograms) and 
simulated (dashed lines) distributions of the radio luminosity at 400 MHz and of distance, 
respectively. We can see that the simulated distribution (dashed lines) qualitatively explain 
the observed features; for example, the peak positions of the distributions. On the other hand, 
we find that it is difficult to explain the sharp peak at 1-2 kpc in the observed distribution 
of the distance. We argue that this excess in the observed distribution is caused by several 
MSPs, and that the excess may be produced by enhancements of the local medium (e.g. the 
Gould Belt) and of the local birth rate. For example, the age of the Gould Belt is believed 
to be ~ 10^^*^ yr (Grenier 2004). Because the birth rate of the canonical pulsars in the 
Gould Belt is expected to be one per 10^ yr, and because the birth rates of the MSPs in the 
Galactic field are 2-3 orders smaller than that of the canonical pulsars, one may expect one 
per 10''"* yr as the birth rate of the MSPs in the Gould Belt. This indicates that several 
MSPs can be born within the age of the Gould Belt. If one assumes that the typical radial 
velocity of MSPs is (70/-\/3) km/s ~ 40 km/s, then the displacement from the Gould Belt 
is (i ~ 5 X 10^ yr x 40 km/s ~ 2 kpc, which is consistent with the enhancement of the peak 
in the distribution of the distance. 

Figures [3(c) and (d) show the averaged luminosity and spin down power, respectively, 
as a function of the distance. The solid and dashed lines are the results for the observations 
and for the simulations, respectively, and each bin of the histograms has the same number 
of observed MSPs. We find in Figures [3t^c) and (d) that the present simulation describes 
characteristics of the observations very well, that is, the averaged radio luminosity increases 
with the distance and the averaged spin down power is almost constant with the distance. 
Finally, we can see that the simulated MSPs can reproduce the Galactic longitude and 
latitude distributions of the radio MSPs. On these grounds, we conclude that our simulation 
reproduces the Galactic distribution of the radio MSPs. 

3.2 Comparison with Fermi observations 

In this section, we compare the results of our simulation with the Fermi six-month long 
observations. Fermi has found 9 "radio-selected" 7-ray emitting MSPs with a fiux F^ ^ 
10~^^ erg/cm^s (Abdo et al 2009b; Abdo et al. 2010a, c), and all of them are located at high 
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Galactic latitudes, that is, \h\ ^ 5°. Abdo et al. (2010a) show that the sensitivity of the low 
Galactic latitudes is about a factor of three worse as compared with that of high Galactic 
latitudes. In the present simulation, therefore, we use F = 10^^^ erg/cm^s for the Galactic 
latitudes |6| ^ 5° and 3 x 10~^^ erg/cm^s for \b\ ^ 5° as the sensitivity requirement of the 
radio-selected MSPs for the six-month long observations. Because no 7-ray-selected MSPs 
have been detected so far, we cannot simulate the Fermi sensitivity of the blind search for 
the MSPs. In this section, we show the results of our simulations by setting the sensitivity to 
that of the canonical pulsars, that is, F = 2 x 10~^^ erg/cm^s for Galactic latitudes \b\ ^ 5° 
and 6 x 10~^^ erg/cm^s for \b\ ^ 5°. 

Table 1 summarizes the population of the radio-selected and the 7-ray-selected MSPs 
detected within the simulations. A^^ and Ng represent the population of the radio-selected 
and the 7-ray selected MSPs, respectively, detected with the Fermi six-month sensitivity 
and Ni^FyiQ^w {i = r,g) are the populations of "bright" MSPs with a 7-ray flux larger than 
IQ-io erg/cm^s. In addition, s and A^i,io are the populations with sensitivity projected 
to 5 year and 10 year observations, respectively, for which the sensitivity is scaled oc v^, 
where T is the length of the observation time. 

As Table 1 shows, our model predicts that only 2-3 radio-selected and 1-3 7-ray-selected 
MSPs can be detected with a flux larger 10"^'' erg/cm^s. Although the results of the Fermi 
with the six-month long observations (Abdo et al. 2009a) did not include the bright sources, 
it is possible that with the small population of the bright MSPs, Fermi may have missed 
the detection of the pulsation. The Fermi unidentified catalog includes about 50 sources 
with a flux larger than 10~^° erg/cm^s. We can see that the spectral behaviors of most of 
them are consistent with those of the known 7-ray pulsars; (1) they are steady sources in 
the sense that the variability index defined in Abdo et al. (2010b) is smaller than 23, and 
(2) the spectra above 100 MeV cannot be fit well by single power law, that is, the curvature 
index, C (defined in Abdo et al. 2010b), is ^ 10. Although most of them may be canonical 
7-ray pulsars (c.f. Takata et al. 2011) or Active Galactic Nuclei (Abdo 2010c), it is very 
likely that a small fraction, in particular, high Galactic sources, are 7-ray emitting MSPs 
(see section |3^ . Recently, Ransom et al. (2011) reported the discovery of three radio MSPs 
associated with the Fermi unidentified bright sources, suggesting that our results are quite 
consistent with the observations. 

Applying the sensitivity of Fermi six-month long observation, our model predicts 9-13 
radio-selected 7-ray MSPs, which is consistent with the 9 MSPs with Fermi observations. 
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We also predict 22-35 7-ray-selected MSPs, whereas Fermi has had no such detections so 
far. For the canonical pulsars, the sensitivity of the blind search is about a factor of 2 worse 
as compared with that of the radio-selected pulsars (Abdo et al. 2010a). The population 
of the 7-ray-selected pulsars is roughly comparable with that of the radio-selected 7-ray 
pulsars. On the other hand, it is expected that the detection of the rotation period from the 
MSPs (in particular in binary systems) by the blind search is even harder and the sensitivity 
is much worse than that of the canonical pulsars. 

Figure H] compares the cumulative distributions of the various characteristics of the 
observed 9 (solid lines) 7-ray MSPs detected by Fermi (Abdo etal. 2009b; Abdo et al. 
2010a,b,c) and the simulated (dashed lines) radio-selected 7-ray MSPs. The results are for 
fmax = fcrit- We performed a Kolmogorov-Smirnov (KS) test to compare the two cumu- 
lative distributions. In Figure HI we present the maximum deviation (Dmax) between the 
two distributions and the p-value (P^s) of the KS-test for each pulsar characteristics. Since 
the simulated sample has more than 400 pulsars, we used the one-sample KS statistic. For 
example, the p-value for the rotation period is Pks ~ 0.75, suggesting that the hypothesis 
that the two distributions are drawn from the same distribution cannot be rejected at better 
than 25% confidence level. We find in Figure H] that for all distributions except that for the 
flux distribution, the hypothesis cannot be rejected at better than 60%, indicating that the 
model distributions are consistent with the observations. 

For the 7-ray flux (right-bottom panel in Figure Hj), the result of the KS-statistic provides 
a relatively low p-value, as compared with other characteristics. This is because all 9 MSPs 
detected by Fermi have a 7-ray flux smaller than 10~^° erg/cm^s, while the model predicts 
that about 20% of 7-ray MSPs have a flux larger than 10~^° erg/cm^s. However, we note that 
the difference between the observed and simulated distributions is caused by only one or two 
MSPs. Hence, if several MSPs with a bright 7-ray emissions of F^ ^ 10^^*^ erg/cm^s have 
been missed by the Fermi observations, the model distribution is not in conflict with the 
observations. It is important to note that our simulation results are based on the existing 
radio sensitivities of pulsar surveys. Any deep radio search implies increasing the radio 
sensitivities and consequently some gamma-ray selected MSPs can become radio-loud. In 
fact, several radio MSPs have been detected by deep search for the radio emissions from 
bright Fermi unidentified sources (Ransom et al. 2011). These radio- loud gamma-ray MSPs 
are defined as gamma-ray selected MSPs in our simulations. 

As Table 1 shows, the model predicts that the population of radio-selected 7-ray MSPs 
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increases by only about 10 over the 5-years (or 10-years) of Fermi observations. This imphes 
that most of the presently known radio MSPs (~ 80) might not be discovered by Fermi. 
However, the predicted radio-loud 7-ray MSPs depend on the sensitivity of the radio surveys, 
as we can see in Table 2. For the 7-ray selected MSPs, the simulation predicts that Fermi 
can detect at least 100 sources. As we will discuss in section 13.41 we expect that these 
simulated 7-ray selected MSPs correspond to the Fermi unidentified sources. As we have 
emphasized before, we predict that more and more radio MSPs will be confirmed among the 
Fermi unidentified sources by a deep search of the radio emissions. 

Story et al. (2007) studied the population of 7-ray MSPs by using the slot gap accelerator 
model and predicted that Fermi would detect about 12 radio-loud and 33-40 radio-quiet 7- 
ray MSPs. On the other hand, our outer gap model tends to predict more 7-ray (in particular 
radio-quiet) MSPs than the slot gap model. Also, the present simulation predicts a larger 
ratio between the radio-quiet and radio-loud MSPs (~ 6 — 7) than that predicted by the slot 
gap model (~ 3), although the ratio depends on the simulated sensitivities of the Fermi 
and of the radio observations (sections 3.3 and 3.4). 

3.3 Population of 7-ray MSPs 

In Figure El we summarize the population of the 7-ray emitting MSPs detected within the 
simulation as a function of the threshold energy flux of the 7-rays. The solid lines and 
dashed lines represent the population of the radio-selected and 7-ray-selected 7-ray MSPs 
respectively. The thick and thin lines represent results for the maximum fractional gap 
thickness of fmax = fcrit and fmax = 1, respectively. For example, the present simulation 
predicts that ~ 2 — 3 and ~ 10 — 14 radio-selected 7-ray MSPs irradiate the Earth with 
a 7-ray flux Fj ^ 10^^° and F^ ^ 10~^^ erg/cm^s, respectively. For the radio-quiet 7-ray 
MSPs, about 100 sources with a flux larger than 10~^^ erg/cm^s irradiate the Earth 

We can see in Figure [5] that the numbers of 7-ray-selected MSPs increase more rapidly 
than that of the radio-selected 7-ray MSPs and the ratio of the 7-ray-selected to radio- 
selected 7-ray MSPs increase with decreasing the threshold energy flux. As the threshold of 
7-ray energy flux decreases, the 7-ray emissions from more distant MSPs can be detected 
in the simulations. On the other hand, because we fix the sensitivity of the radio surveys 
in Figure O the radio emissions from those distant MSPs may not be detected by the 
radio surveys. As a result, more 7-ray selected MSPs are detected in the simulations as the 
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threshold flux decreases. In fact, if we count all 7-ray MSPs irradiating the Earth with the 
radio emissions as the radio-selected pulsars, the ratio does not depend on the threshold 
energy flux of the 7-ray emissions. 

Figure [6] shows contours of the two-dimensional distribution in the rotation period and 
the period-time derivative of the simulated MSPs with a 7-ray flux ^ 10^^^ erg/cm^s. 
We present the distributions for all simulated MSPs including the radio-selected and 7-ray- 
selected pulsars, because the two distributions do not differ much from each other. The left 
panel shows the distributions using the period time derivative for the intrinsic value, and the 
right panel represents the period time derivative for the "observed" value, that is, the value 
after adding the Shklovskii effect described by equation ( IT2l) . Within solid, dashed and dotted 
lines, 10%, 50% and 100% of the total MSPs are populated. We can see that the Shklovskii 
effect slightly shifts the peak of the distributions toward longer rotation periods and larger 
period time derivatives. For the observational distributions (right panel), about 50% of the 
7-ray emitting MSPs will be detected with a rotation period between P ~ 0.003 — 0.008 s 
and a period time derivative between P ~ 6 x 10~^^ — 10"^^. This information can be used 
to narrow down the parameter range for a blind search to detect the rotation period from 
Fermi unidentified sources. 

3.4 Origin of the Galactic Fermi unidentified steady sources 

As Figure |5] shows, our simulation predicts that about ~ 200 7-ray MSPs irradiate the Earth 
with a flux F^ ^ 5 X 10~^^ erg/cm^s and with the present sensitivity of the radio surveys 
most of them are categorized as 7-ray-selected MSPs. Therefore, it is likely that although 
millisecond pulsed emission has not been confirmed yet, the 7-ray MSPs contribute to the 
galactic unidentified Fermi sources, such as the newly discovered 20 radio MSPs associated 
with the Fermi unidentified sources (Ray 2010; Caraveo 2010; Ransom et al. 2011). 

Figure [7]plots the curvature index (C) and the variability index (V) of the Fermi sources; 
a curvature index larger than 11.34 indicates a less than 1% chance that the power-law 
spectrum is a good fit above 100 MeV, and a variability index larger than 23.21 indicates 
less than a 1% chance of being a steady source (see Abdo et al. 2010b for the exact definition 
of C and V indexes). The filled circles, boxes and triangles correspond to the pulsars. Active 
Galactic Nuclei (AGN) and Fermi unidentified sources. In Figure [71 we can easily see that 
the pulsars and the AGN belong to different groups in the C-V plane, that is, the pulsars 
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belong to C ^ 5 and V ^ 20, whereas AGN belong to wide range of variability index, which 
may indicate a correlation between the C- index and V-index (that is, C ocV). Most of the 
Fermi unidentified sources have the V-index smaller than ~30, with the C-index larger than 
0.1. Figure [7] indicates that the Fermi unidentified sources can be either pulsars or AGNs. 

In Figure [H we plot the Galactic distributions of the Fermi unidentified steady sources 
for C < 5 (solid line) and C ^ 5 (dashed line). It is clear from Figure [8] that the two lines 
represent different distributions with respect to each other. The solid lines show relatively 
constant distributions for the Galactic longitudes and latitudes, indicating that most of the 
Fermi unidentified sources with C < 5 (dashed lines) distribute isotropically in the sky, 
and may be related with extra-Galactic (or Galactic halo) sources. For the dashed-lines, 
on the other hand, the distributions for the Galactic longitudes and latitudes coordinate 
show peaks at the Galactic centre and at the Galactic disk, respectively, suggesting that the 
Fermi unidentified steady sources with C ^ 5 are associated with Galactic objects. 

In Figure [HI we compare the distributions of the Galactic longitudes (left panel) and 
latitudes (right panel) for the unidentified Fermi sources with V ^ 23.21 and C ^ 5 (solid 
line) and the simulated 7-ray MSPs (dashed line), with a flux ^ 10"^^ erg/cm^s. The 
latitude distribution of the 7-ray emitting canonical pulsars (dotted line) simulated in Takata 
et al. (2011) and the Galactic distributions for known radio MSPs (dotted lines) are also 
plotted in the flgure. 

In Figures [9] we can see that the distributions of the simulated 7-ray MSPs and of the 
observed radio MSPs are consistent with that of the Fermi unidentified sources, that is, (1) 
three longitude distributions have a peak around the direction of the Galactic centre / = 0° 
and become minimum around / ~ 180°, and (2) the latitude distributions have a peak at 
the Galactic plane (6 = 0) and then they decreas with the increase of the Galactic latitudes. 
In particular, the MSPs can explain the distributions of the Fermi unidentified sources 
located above the Galactic plane \h\ ^ 10, which cannot be explained by the canonical 7-ray 
pulsars, as the right panel in Figure |9] shows; the 7-ray emitting canonical pulsars can mainly 
explain the unidentified sources located in the Galactic plane. Since the MSPs are in general 
older than the canonical pulsars, a higher fraction of the 7-ray MSPs, as compared with the 
canonical pulsars, is located at higher Galactic latitudes. On these ground, we conclude that 
7-ray emitting MSPs are more plausible as candidates for the origin of the majority of the 
Galactic Fermi unidentified steady sources located in high Galactic latitudes. 

We note that new radio MSPs have been discovered in the direction of about 20 uniden- 
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tified Fermi sources by the radio searches for individual source (e.g. Ray 2010; Caraveo 
2010; Ransom et ah 2011). These '-/-ray- selected radio-loud MSPs can be classified as the 
7-ray-selected radio-quiet sources in the present simulation, because we do not take into 
account the radio search for the specific sources. In Table 2, on the other hand, we present 
how the population of 7-ray emitting MSPs depends on the sensitivity of the radio survey. 
We used the Arecibo 2 (A2) and Parks 2 (P2) surveys (the first row), all radio surveys (the 
second row) listed in table 1 of Takata et al. (2011) and all radio survey but we increase by 
a factor of 2 the sensitivity of each survey (third row). The intrinsic populations, that is, 
the populations associated with only beaming effects of the radio emission are shown in the 
bottom row. 

We can see in Table 2 that the number of radio-selected pulsars increases with the 
increase of the sensitivity of the radio surveys, whereas the 7-ray-selected pulsars decrease. 
This is because if we increase the sensitivity of the radio surveys some 7-ray-selected MSPs 
irradiating the Earth with the radio emissions may be re-classified as radio-selected 7-ray 
pulsars. As the bottom row in Table 2 shows, the intrinsic ratio of the radio-loud and radio- 
quiet 7-ray MSPs is 10~ 20, indicating that most of the 7-ray emitting MSPs irradiate 
the Earth with a corresponding radio emission. Therefore, our prediction is that as the 
sensitivity of the radio observations will improve in the future, more and more radio MSPs 
will be discovered among the Fermi unidentified sources. 

In Figure [101 we present the distributions of various characteristics of observed 7-ray 
MSPs (shaded histograms), including the 20 new radio MSPs associated with Fermi uniden- 
tified sources, and the simulated 7-ray MSPs (dashed lines) with the simulated sensitivity 
of the Fermi six-month long observation. For the simulated distribution, both the radio- 
selected and 7-ray-selected MSPs are taken into account. For the observations, we were able 
to obtain information of the rotation period for 7 out the 20 radio MSPs and the distances 
for 3 MSPs (Ray, 2010; Ransom et al. 2011). Therefore the distributions of the rotation 
period and of the distance in Figure [TOl are obtained from 16 (9 known Fermi MSPs plus 7 
radio MSPs) and 12 samples, respectively. The distributions of the period time derivative, 
the magnetic field and the spin down age include only 9 known Fermi MSPs, because we 
could not find any references for the intrinsic period time derivative of the 20 radio MSPs. 
The 7-ray fiuxes of the 20 radio MSPs are taken from the Fermi first catalog (Abdo et al. 
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2010b 0) and Ransom et al. (2011). In Figure fTOl we also indicate the p-values {Pks) of the 
KS test. As Figure ITO] shows, the simulated distributions are qualitatively consistent with the 
observations (e.g. the position of the peak of the distributions). For the 7-ray flux, however, 
the p-value of KS-test is much lower than those of other characteristics, indicating that the 
two distributions cannot be drawn from same distribution. Because it is expected that the 
7-ray emissions from the 20 radio MSPs associated with Fermi unidentified sources will be 
composed of a pulsed and an unpulsed component, the extraction of the pulsed component, 
that is, the detection of the pulsed period with Fermi, is required to constrain our model. 

We have assumed that MSPs are activated as the rotation powered pulsar as a result of 
the accretion process from the companion star. Recently Takata, Cheng & Taam (2010b) 
argued that the outer gap activities of MSPs can be turned on in the quiescent state of low 
mass binary systems, for example PSR J1023-I-0038 (Archibald et al. 2009). The geometry 
of the outer gap can direct the gamma-iajs emitted from the gap towards the companion 
star and the accretion disk. Consequently, the accretion disk can be evaporated, and the 
companion star becomes brighter due to the irradiation of gamma-iajs. If the high latitude 
unidentified sources of Fermi are indeed MSPs, some of them may be low mass binaries in 
a quiescent state. We expect them to be associated with abnormally bright stars with short 
orbital period. By searching the orbital periods of the optical companion stars, we may be 
able to obtain the spin periods of the radio-quiet MSPs from gamma-ray data. 

4 CONCLUSION 

We have studied the population of the 7-ray emitting millisecond pulsars using Monte- 
Carlo techniques. We have applied the outer gap model with a switching of the gap closure 
process from the photon-photon pair-creation model to the magnetic pair-creation model, 
as suggested by Takata et al. (2010a). Using the sensitivity of the Fermi six-month long 
observations, 9-13 radio-selected 7-ray MSPs are detected within the simulation, which is 
consistent with the present 9 Fermi MSPs. The simulated distributions for the various 
characteristic of the radio-selected 7-ray pulsars are consistent with the present Fermi 
observations (Figure H]). Scaling the observed sensitivity oc \/T, the present model predicts 
that Fermi should detect 15-22 radio-selected 7-ray MSPs and 95-152 7-ray-selected MSPs 
over its 5- year mission (Table 1). 



see also http : //fermi. gsfc.nasa.gov/ssc/data/access/lat/lyr_catalog 
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Our simulation predicts that about 100 (or 200-300) 7-ray MSPs irradiate the Earth 
with a flux F ^ 10"^^ erg/cm^s (or 5 x 10^^^ erg/cm^s), and most of them are categorized 
as 7-ray-selected pulsars with the present sensitivity of the radio surveys (Figure [5]). Our 
simulation also predicts that about 50 % of the 7-ray emitting MSPs will be detected with 
a rotation period in the range P ~ 3 — 8 ms and a period time derivative in the range 
P ~ 6 X 10^^^ — 10^^^. We further argue that 7-ray emitting MSPs are plausible as candidates 
for the Galactic Fermi unidentified steady sources, located in high Galactic latitudes. Our 
simulation implies that some of the radio-quiet 7-ray MSPs can be changed to radio-loud 7- 
ray MSPs as the simulated sensitivity of the radio surveys increasees. Therefore, we predict 
that more and more radio MSPs will be discovered in the Fermi unidentified sources as the 
radio sensitivity is improved. 
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Table 1. Population of the simulated radio-selected and 7-ray-selected 7-ray MSPs. The subscript r and g represent the 
radio-selected and 7-ray-selected pulsars, respectively. pyio-^o ~ d) ^^'^ population of "bright" MSPs with a 7-ray flux 
larger than 10"^" erg/cm^s, and A^^ are results with the Fermi six-month long observation. In addition, Ni^^ and Ni^io are 
population with the sensitivity projected to 5 yr and 10 yr observations, respectively. 
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Table 2. Population of the simulated radio-selected and 7-ray-selected MSPs for the different radio surveys and projected 
mission lengths of Fermi(6 months, 5 years and 10 years). A2 and P2 represent the Arecibo 2 and Parks 2 surveys. "AU" 
represents all radio surveys listed in table 1 of Takata et al. (2010b) and "Allx2" shows the results with all radio survey but we 
increase a factor of 2 for the sensitivity of each survey. In addition, the intrinsic populations, that is, the populations associated 
with only beaming effects of the radio emission, are shown in the bottom line. In this extreme limit, most Fermi 7-ray pulsars 
become radio-selected 7-ray pulsars. 
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Figure 1. The observed distributions (shaded histograms) of the various characteristics for radio MSPs. The dashed linos 
show the distributions applied for the present simulation. The observed distributions are for the intrinsic values (i.e. without 
Shklovskii effect) of the pulsar parameters. 
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Figure 2. The fractional 7-ray flux as a function of tfic viewing angle and of the inclination angle. The fractional flux is defined 
by the ratio of the calculated flux measured on the Earth to f^Lsdl^^ ■ The results are for P = 3 ms, _Bs = 3 X 10® Gauss and 
f(Ri,) = 0.4. 
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Figure 3. Distributions of (a) luminosity and (b) distance of the observed (shaded histograms) and simulated (dashed lines) 
radio MSPs. The panels (c) and (d) represent the averaged radio luminosity and the spin down power for the observed (solid 
lines) and the simulated (dashed lines) MSPs as a function of the distance, respectively. Each bin of the solid lines contains 
the same number of the observed MSPs. 
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Figure 4. Cumulative distributions for the various characteristics of the observed (sohd Uncs) and simulated (dashed lined) 
radio-selected 7-ray MSPs with the sensitivity of the Fermi six- month long observation. The maximum difference (Dmax) and 
the p- value of KS-test (Pks) E're also displayed. The intrinsic values, that is, the values without the Shklovskii effect, are also 
represented. 
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Figure 5. The simulated numbers of the radio-selected (solid line) and 7-ray-selected (dashed line) 7-ray MSPs as a function 
of the threshold energy flux. 
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Figure 6. Contours of two-dimensional (rotation period and period time derivative) distributions of the simulated MSPs with 
the sensitivity of the Fermi six-month long observation; (a) intrinsic period time derivative and (b) "observed" period time 
derivative adding Shklovskii effect. More MSPs are located within the solid lines. Within solid, dashed and dotted lines, 10%, 
50% and 100% of the total MSPs are populated. 
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Figure 7. Plots for the curvature index and tlie variability index of the Fermi sources (Abdo et al., 2010b). The filled circles, 
and boxes and triangles are the pulsars, AGNs, and the Fermi unidentified sources respectively. 
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Figure 8. The distributions in Galactic coordinate longitudes (left) and latitudes (right) of the Fermi unidentified steady 
sources. The solid and dashed lines arc the results for C < 5 and C ^ 5, respectively. 
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Figure 9. The distributions of Galactic longitudes (left) and latitudes (right) of the Fermi unidentified point sources with 
V ^ 23.21 and C ^ 5 (solid lines), the simulated 7-ray MSPs (dashed lines) and observed radio MSPs (dotted lines). The 
distribution of Galactic latitudes of the simulated 7-ray emitting canonical pulsars (CNP) in Takata ct al. (2010b) are also 
plotted in the right panel (thin dotted line). 
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Figure 10. Distributions of the various properties of the observed (shaded histograms) and of the simulated (dashed lines) 
MSPs. In addition to the 9 known 7-ray MSPs, the observed distributions of the rotation period and of the distance includes 7 
and 3 radio MSPs, respectively, associated with the Fermi unidentified sources. For the 7-ray flux, the observed distributions 
includes 20 radio MSPs associated with the Fermi unidentified sources. For the simulation, both radio- selected and 7-ray- 
selectcd MSPs are included in the distributions. 



